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Abstract. We consider two protoplanets gravitationally 
interacting with each other and a protoplanetary disc. The 
two planets orbit interior to a tidally maintained disc cav- 
ity while the disc interaction induces inward migration. 
When the migration is slow enough, the more rapidly mi- 
grating outer protoplanet approaches and becomes locked 
in a 2 : 1 commensurability with the inner one. This is 
maintained in subsequent evolution. We study this evolu- 
tion using a simple analytic model, full hydrodynamic 2D 
simulations of the disc planet system and longer time N 
body integrations incorporating simple prescriptions for 
the effects of the disc on the planet orbits. The eccentric- 
ities of the protoplanets are found to be determined by 
the migration rate and circularization rate induced in the 
outer planet orbit by the external disc. 

We apply our results to the recently discovered reso- 
nant planets around GJ876. Simulation shows that a disc 
with parameters expected for protoplanetary discs causes 
trapping in the 2 : 1 commensurability when the plan- 
ets orbit in an inner cavity and that eccentricities in the 
observed range may be obtained. 
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1. Introduction 

The recent discovery of extrasolar giant planets orbiting 

around nearby solar-type stars (Marcy & Butler 1998, 

2000) has stimulated renewed interest in the theory of 

planet formation. The objects observed so far have masses, 

Mp, that are characteristic of giant planets 

(0.4 Mj < Mp < 11 Mj), Mj denoting a Jupiter mass. 

The orbital semi-major axes are in the range 

0.04 AU ^ a < 2.5 AU, and orbital eccentricities in the 

range 0.0 < e < 0.67 (Marcy & Butler 2000). 

Disc-protoplanet interactions have been invoked to 
explain the presence of giant planets orbiting close to 
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their host stars through inward orbital migration in- 
duced through disc-protoplanet tidal interaction ( eg. Pa- 
paloizou, Terquem & Nelson 1999, Lin et al 2000). Up to 
now, for the most part extrasolar planets appear to be 
isolated. However, a few multiple systems are known. The 
configuration of these may contain important informa- 
tion about their origin and possible migration history. Of 
special interest is the recently discovered system around 
GJ876. This is found to be close to a 2 : 1 commensura- 
bility. Such a configuration is indeed suggestive of orbital 
migration. Commensurable satellite systems such as the 
Galilean satellites are thought to owe their origin to migra- 
tion induced by tidal interaction with the central planet ( 
eg. Goldreich 1965). 

Recent simulations of single protoplanets in the ob- 
served mass range (Kley 1999, Bryden et al. 1999, Lubow, 
Seibert & Artymowicz 1999) interacting with a disc with 
parameters thought to be typical of protoplanetary discs, 
but constrained to be in circular orbit, indicate gap for- 
mation and upper mass limit consistent with the observa- 
tions. Nelson et al. (2000) relaxed the assumption of fixed 
circular orbits, found inward migration and that the disc- 
protoplanet interaction leads to strong eccentricity damp- 
ing. Due to accretion onto the central star, an inner cavity 
was formed in the disc interior to which the protoplanet 
orbited. 

Simulations of two planets interacting with a disc have 
been performed by Kley (2000), Bryden et al (2000), and 
Masset & Snellgrove (2001). So far inward migration of 
two planets locked into a 2 : 1 commensurability has not 
been simulated. However, Bryden et al (2000) found a ten- 
dency for gap material between the two planets in fixed 
circular orbits to be cleared, ending up interior to the in- 
ner planet orbit or exterior to the outer planet orbit. 

Taken together the above results suggest a natural out- 
come of two protoplanets interacting with a disc is that 
they orbit interior to an inner disc cavity while the exter- 
nal disc causes inward migration of the outer orbit. This 
catches up the inner orbit leading to the possibilty of res- 
onant interaction. 

It is the purpose of this paper to investigate such reso- 
nant interaction and whether, for reasonable protoplane- 
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tary disc models, it leads to a locking such that the planets 
subsequently migrate maintaining the commensurability. 
Similar behaviour occurs as a result of the tidally induced 
migration of the Gallilean satellites (eg. Goldreich 1965, 
Lin & Papaloizou 1979). 

Here we shall assume the prior evolution of the system 
leads to the orbital separation of the planets being slightly 
larger than that required for a strict 2 : 1 commensurabil- 
ity without considering the history in detail as it is beyond 
the scope of this paper. However, we comment that this 
might have been complicated with the planet masses vary- 
ing with time through mass accretion from the disc. 

Analytic methods, N body integrations and direct two 
dimensional numerical simulations are used to investigate 
the evolution and found to give consistent results. 

We consider two protoplanets gravitationally interact- 
ing with each other and a protoplanetary disc. The two 
planets orbit interior to a tidally maintained disc cavity 
while the disc interaction induces inward migration. This 
migration is of type H and so is regulated by the mag- 
nitude of the disc viscosity. When the migration is slow 
enough, the more rapidly migrating outer protoplanet ap- 
proaches and becomes locked in a 2 : 1 mean motion com- 
mensurability with the inner one. The commensurability 
persists in subsequent evolution. The eccentricities of the 
protoplanets are increased by the resonant perturbations 
and damped by circularization which occurs through in- 
teraction with the disc (eg. Goldrech & Tremaine 1981). A 
balance is achieved in which the protoplanet eccentricities 
are determined by the migration rate and circularization 
rate induced in the outer planet orbit by the external disc. 

We apply our results to the recently discovered res- 
onant planets around GJ876. Simulation shows that mi- 
gration induced by a disc with parameters expected for 
protoplanetary discs results in trapping in the 2 : 1 com- 
mensurability when the planets orbit in an inner cavity. 
Eccentricities in the observed range may be obtained. Fur- 
ther studies using N body integrations indicate that the 
planetary system will remain stable for at least 2 x 10^ 
orbits when the external disc is removed. 

In section we describe a simple analytic model of 
two migrating protoplanets in a 2 : 1 commensurability. 
The eccentricities of the protoplanet orbits are related to 
the migration rate, and circularization rate induced by the 
disc. In section (||) we describe a simulation of two planets 
orbiting in an inner disc cavity. Parameters appropriate to 
GJ867 are adopted. This demonstrates resonant trapping 
and that eccentricities of the observed magnitude may be 
produced. 

In section (Q) we describe N body calculations confirm- 
ing the above conclusions and indicating the long term 
stability of the system. Finally in section (|^) we discuss 
our results. 



2. A simple model 

We consider a system consisting of 2 planets and a primary 
star moving under their gravitational attraction. When 
there are no disc interactions and the motion is conserva- 
tive the system is conveniently expressed in Hamiltonian 
form using Jacobi coordinates (eg. Sinclair 1975). The co- 
ordinates, r2, of inner planet of reduced mass rn2 are re- 
ferred to the central star of mass Af, and the coordinates 
of the outer planet, ri, of reduced mass mi are referred 
to the centre of mass of the central star and inner planet. 
The Hamiltonian can be written correct to second order 
in the planetary masses as 



H 



1 
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Here A/*i — Ah + mi , M*2 = Ah + m-i and ri2 = r2 — ri . 
The Hamiltonian can be expressed in terms of the os- 
culating semi-major axes, eccentricities and longtitudes 
of periastron ai,ei,Wi,i — 1,2 respectively as well as 
the longtitudes Aj, and the time t. We recall that Aj = 
^iit — tgi) + Wi, with rii being the mean motion and tgi 
giving the time of periastron passage. The energy is given 
by Ei = —G'miA'hi/(2ai), and the angular momentum 
hi — ■mi^GAhiai{\ — ef) which may be used to describe 
the motion instead of at and e^. 

Only terms first order in the eccentricities and in- 
volving the resonant angles 4> = 2Ai — A2 — wi, and 
tp — 2Ai — A2 — 1^2, are retained in the expansion (see 
Sinclair 1975). Then the perturbing part of the Hamilto- 
nian (oc mim2) can be written 



= {Bei cos(f> — Ce2 cosip) 



ai 



(2) 



where C = 2b' 



[J^{a) + {l/2)adb[y^{a)/da and B 



(2) 



1. 



.56^/1,(0!) + {l/2)adbY/2{a)/da 
the usual Laplace coefficient and a — 02/01. From now on 
we replace A'hi by M*. 



(1) 



2a. Here 6^"2 denotes 



2.L Orbital precession 

We may also take into account additional gravitational 
forces that may produce precession of the planetary orbits. 
These could result from the mass distribution provided by 
the disc exterior to the outer planet or the time averaged 
mass distribution of the planets themselves. These effects 
are not included in the resonant Hamiltonian given above. 
Accordingly we add to it another Hamiltonian 



Hint = ~-jminia{LOprie\ - ^m2n2a\u}pr2e\- 



(3) 



Here, as can be verified from the equations of motion, we 
have adopted a parameterization such that the precession 



Snellgrove, Papaloizou & Nelson: Resonance coupling of planets in a disc 



3 



frequency induced in the orbit of rrii is Upri^ai). This pre- 
scription does not allow the precession frequecy to depend 
on other quantities. However, appropriate matching can 
be carried out for a particular case under consideration. 

2.2. Basic Equations 

The equations of motion are derived from: 
dEi/dt ^ -^UidH'/dXi - (niT + D)Sii, 
dhjdt = -dH'/dX, - dH'/dwi - TSu, 
d\,/dt = n, + UidH'/dE, + dH'/dhi, 
dwi/dt = dH' /dhi, 

with H' = Hi2 + Hint. These can be obtained from Hamil- 
ton's equations (eg. Brouwer & Clemence 1961) to which 
we have added, for the outer planet mi, an additional ex- 
ternal torque —T with an associated orbital energy loss 
rate niT together with additional orbital energy dissipa- 
tion rate D. The torque and dissipation rate could be pro- 
duced by tidal interaction with the disc leading to inward 
migration and orbital circularization. 

We thus obtain to lowest order in the planetary eccen- 
tricities and perturbing masses. 



dni 

~dr 



dn2 

~dr 

dei 
~dt 

de2 
~dt 



GM^mi 
?tn\mia2 



{Bei sin 4> — Ce2 sin ip) 

{niT- 



D) 



{Bei sin0 — Ce2 sin^') 



TO2711 



B sin (j) — 



Da, 



771171202 

M*ai 



C sin-i/) 



7712 



= 2711 ~ ^2 

dt eiAf, 



niB cos (/) — ujpri 



dip _ , 771102 
— = 2711 - 712 + — 

dt e2aiM* 



2.3. Stationary solutions 



7l2CcOS^ — LUpr2 



(4) 
(5) 
(6) 
(7) 
(8) 
(9) 



When no migration or circularization occurs {T — D — 0) 
equilibrium solutions may exist such that ip and (f) are 
either zero or tt. Each of 7ii, 712, ei, 62 are then constant. 

A relation between the eccentricities then follows from 

and (||) in the form 



e2ai{eiM^LUpri + m2niB coscf)) = 
ei (—771102712(7 coai/j + e2aiAf,Wpr2)- 



(10) 



This condition matches the precession rates of the orbits of 
the two planets. Also 27ii = 712. Noting that the eccentric- 
ities are positive, when they are of very small magnitude, 
the precessional terms become negligible and there is a 



solution with -0 
we have 



0, (/) = TT or -0 — 0, (/) = vr. In either case 



7ii2aie27ii_B = 77iia2ei7i2C. 



(11) 



For larger eccentricities the precessional terms may be- 
come important in ( |lO|) a nd then solutions with ip = 0, (j> = 
0, may occur. Then (|lO|) gives 



77l2aie27li_B -I- 771i02ei7l2C — eie20iAf,(cJpr2 



Op^l 



)• (12) 



For stable solutions, when perturbed, the angles may 
undergo librations about their equilibrium points (eg. 
Sinclair 1975). There are two frequencies of oscillation 
1/1,^2 being given for any ip,(j) in the limit of small 
eccentricities by i^f = (m2?^iS)^/(MH,ei)^, and f| = 
(7iii7i2a2C)^/(aiM,e2)^. We look for a solutions with mi- 
gration which are close to stable solutions of this type. 

2.4. Resonant Migration 

We look for solutions of (Qj|) corresponding to the situa- 
tion where the two planets migrate inwards locked in res- 
onance with 7ii/n2 maintained nearly equal to 1/2 while 
the eccentricities remain nearly constant. The tendency 
for the resonant coupling to excite the eccentricities is 
counterbalanced by circularization through the action of 
D = {GM^,mie1) / (aitc) which defines the circularization 
time for ei. Similarly (^ defines an inward migration 
timescale tmig — GAI^mi/ {STaiUi). 

We begin by supposing that the angle tp executes a 
libration about zero such that the mean rate of change 
of 62 is zero. Similarly the mean rates of change of 711 
and 7i2 induced by -0 are zero. Such a libration is seen in 
simulations. We also suppose the angle either librates 
or circulates but in such a way that the correspondingly 
induced mean rates of change are not zero. The simplest 
example is when the angle executes a very small or even 
zero amplitude libration about a value slightly offset from 
zero or tt (eg. Lin & Papaloizou 1979). 

We suppose the circulation/libration periods to be 
short compared to the timescale of migration so that av- 
eraging is possible. We denote the average of 7iiei sin ip by 
5. Then averaging (0^) gives 



1 dni 
ni dt 

1 dn2 
71 2 dt 



67712-6(5 



3oi 



GAf,77ll 



(niT + D) 



dt 



dt 



6771l_B02(S 
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27712-8(5 2I?oi 



Af, 



GM^rui 



1712 



niB cos ( 



7711O2 
e2Af»0] 
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(13) 
(14) 

(15) 
.(16) 



Here on the left hand sides the overline denotes the 
time average and for simplicity of notation this has been 
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dropped from the right hand sides. The resonance condi- 
tion also implies — — — . Using this and eliminat- 

^ n2 at ni at ^ 

ing S from the above, we obtain for the rate of increase of 
rii through migration 



f VS time 



1 dni 

Hi dt 



3ai 



{niT + D) 



TOia2 



mia2 + ^1201 



Also (|l|) gives for the eccentricity balance 
tcm2ai 



el 



(17) 



(18) 



3tmig(2TOia2 + m2ai) 

The above determines the eccentricity of the outer 
planet ei as a function of tc and t„iig. For a system with 
raxlra-i = 3, we get ei ^JO^GTtJUm'g . 
For ei in the 0.01 range we need ^ 10 orbits if tmig ~ 
10^ orbits. 

The eccentricity of the inner planet is determined by 
equation (|l^). For small amplitude librations this is given 
by (p^. That would still apply when (/> is circulating pro- 
vided the cosines are time averaged and the mean circu- 
lation rate is small. 
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3. A simulation of two migrating resonantly 
coupled planets 

The protoplanetary disc is numerically simulated using 
an Eulerian 2D hydrodynamic code. The code used is a 
modified version of NIRVANA, which has been described, 
tested and used successfully elsewhere on a similar prob- 
lem involving interacting planets (Masset & Snellgrove 
2001 ) . Incorparated with the hydrodynamic code is a 4th- 
order Runge-Kutta integrator which is used to evolve the 
orbits of the two planets. The gravitational forces calcu- 
lated from the disk model are used in the equations of 
motion of the planets, and the disc itself responds to the 
planetary potential. Hence the system evolves in a self- 
consistent fashion. In order to obtain the long integration 
times needed f or sim ulations of this type the FARGO algo- 
rithm (Masset 2000) is applied. However, tests have shown 
that the results are not affected by this. 

We use a 2D cylindrical (r, (p) grid with 200 radial 
zones distributed uniformly between r — 0.4 and r — 3.47 
in dimensionless units and 300 azimuthal zones. We apply 
outflow conditions at the inner boundary to simulate the 
accretion of disc material onto the central star. 



3.1. Physical setup 

We attempt to simulate the resonance locking of the sys- 
tem GJ876 via tidally induced migration of the planets, 
using plausible values of the disc parameters. The disc 
is assumed to be thin and isothermal, with constant as- 
pect ratio h/r — 0.07, and a constant Shakura & Sun- 
yaev ( |1973| ) a-viscosity prescription with a = 2 x 10"'^ 
is adopted. The two planets are initially in circular or- 
bits coplanar with the disc, at radial locations ri — 1.0 



Fig. 1. The lower plots show the semi major axes a (lower 
left) and eccentricities e (lower right) for both planets. The 
upper plots show the resonance angles (j) and if). The time 
unit is measured in orbit periods of the disc material at 
r = 1. 



and r2 = 0.6. Hence the outer planet is located outside 
the exact 2 : 1 commensurability (r2 = 0.63). The planet 
masses are chosen to correspond to the minim um m ass ra- 
tios obtained from observations (Marcy et al. ^001 ). With 
masses normalised so that stellar mass is Af, ~ 1, this 
corresponds to mi = 6 x 10"'^ and m2 = 1.8 x 10^^. The 
planet masses are fixed as the planets are assumed to be 
no longer accreting material from the disc. 

The disc is prescribed an initial surface density Eo cor- 
responding to what would give a disc mass of 2 x 10^"^ 
within the orbit of the outer planet. However we as- 
sume that both planets are located inside a tidally trun- 
cated cavity located at r < 1.3, with low surface density 



Scavity = O.OISq. This cavity is supposed to have already 
been cleared by the tidal action of the two planets. Be- 
tween 1.3 < r < 1.5 the surface density is prescribed such 
that InS linearly joins to Eg.. 



3.2. Results 

The tidal interaction of the planets with the disc material 
causes the planets to migrate inwards (see Figure |^) . The 
inner planet is deep within the cavity and only interacts 
with low surface density material and thus migrates slowly. 
The outer planet has its outer 2 : 1 lindblad resonance 
located outside the cavity and in the body of the disc. 
Hence there is more material exerting a negative torque 
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Surface density contour plot 




Fig. 2. Surface density plot for inner regions ofthe disc 
at a time corresponding to 300 orbits. Darker areas on 
the plot correspond to regions of lower surface density. 
The disc cavity is clearly visible, as are the density waves 
outside the cavity. Inside the cavity are the wakes of the 
two planets. 



on the planet, and therefore it migrates faster, despite 
its larger mass. The ratio of semi- major axes 01/02 of 
the planets decreases until the planets 'lock' into a 2 : 1 
commensurability with n2 ~ 2ni at a time t « 400 orbits. 
Both planets then subsequently migrate inwards a further 
10% maintaining this ratio, showing the resonance to be 
robust. 

Figure |l] also shows the calculated values of the reso- 
nance angles (j) and ij}- Once the commensurability lock has 
occured, these are both librating about zero. The resonant 
interaction causes a eccentricity growth of both planets, 
the growth halting at around average values of ei = 0.06 
and 62 = 0.34 although both eccentricities exhibit varia- 
tions around these average values. The perihelion angles 
of the two planets oscillate around the alignment position, 
being the natural stable state. 

The disc cavity remains at low density, and the cavity 
edge diffuses slowly in on the viscous diffusion timescale, 
following the outer planet. Figure |^ shows a surface den- 
sity plot of the inner regions of the disc. 

3.3. Relation to analytic model and observations 

The simulation confirms that tidally induced migration 
of two planets with the given mass ratios leads to trap- 



ping into the 2 : 1 commensurability. The resonance angles 
both undergo stable periodic librations as per section ^ 
and the resonance is robust for the length of the simu- 
lation. We can conclude that stable resonance trapping 
would be the probable outcome of the evolution of such 
a system. The finding of the system GJ876 near to or in 
such a resonance is probably simply a consequence of past 
tidally induced migration of the two planets into such a 
state. The magnitudes of the eccentricities are in broad 
agreement with fitted model parameters to the observa- 
tio ns of GJ876 with ei w 0.1 and 62 « 0.27 (Marcy et 
al. 2001). The eccentricity ratio 62/61 ~ 5.7. We com- 
ment that equation (|l^) which applies in the small ec- 
centricity limit when one of 0, ip is zero and the other 
TT gives ei/e2 ~ 1/11. However, the fact that both angles 
librate about zero indicates, within the context of the sim- 
ple model, that non resonant orbital precession needs to 
be incorporated and we should use equation ( p^ to make 
a comparison with the simulation. This is 



616201 A/^Wpri -I- rn20i62ni_B cos — 
— r7ii026in2C COS?/' + 6i620iM*a;pr2- 



(19) 



Since the librations of the angles are about a value close to 
zero, we replace them by zero to make simple estimates. 
The largest non central mass in the system is mi so we 
include its effect in causing non resonant orbital precession 
for m2 by including a non zero ujpr2- Also we include the 
corresponding effect on mi due to m2 which produces a 
non zero value of Wpri • Then we have 



m2niB 712^,102(7 

Wprl H 77 — = I^pr2 — 



61M* 



620iAf, 



(20) 



This condition equates the precession rate of the orbit 
of m2 on the right with that of mi on the left. In the 
former case the first contribution comes from the time av- 
eraged orbit and the second is the resonant contribution. 
In fact with 62 ~ 0.34, the orbits approach each other 
quite closely to within 22 percent of the outer semi major 
axis. But the resonant configuration avoids close encoun- 
ters keeping the planets apart. Thus we expect the reso- 
nant and non resonant contributions to the precession rate 
to show significant cancelation. The large value of 62 prob- 
ably makes strict comparison innacurate. Nonetheless the 
simulation gives a precession period of about 70 orbits in 
a retrograde sense. We can estimate that this to be of the 
same magnitude as would be predicted from the second 
term on the left hand side which however gives prograde 
precession. Thus the non resonant precession rate due to 
m2 needs to be be comparable to the resonant effect but 
of opposite sign. 

Considering the precession rate of m2, we estimate that 
ujpr2 corresponds to 70 orbits in the prograde sense while 
the resonant term corresponds to 70/2 orbits in the retro- 
grade sense. Thus the combined effect produces a preces- 
sion period of 70 orbits in the retrograde sense as seen in 
the simulation. 
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Run 




to 




ei 


62 


62/61 


Rl 


3.33 X 10^ 


900 


8 X 10^ 


0.095 


0.41 


4.3 


R2 


3.33 X 10^ 


250 


8 X 10^ 


0.05 


0.3 


6 


R3 


3.33 X 10^ 


900 


8 X 10^ 


0.34 


0.72 


2.1 


R4 


10^ 


270 


2.4 X 10^ 


0.095 


0.41 


4.3 


R5 


3.33 X 10^ 


900 


0.0 


0.095 


0.41 


4.3 



Eccentricity v's Time 



Table 1. This table shows the parameters used for the 
orbital integration runs. The first column gives the run 
label, the second column the migration time scale, the 
third column the circularisation time scale, and the fourth 
column gives the disc dispersal time. The fifth, sixth and 
seventh columns give the final eccentricities for the outer 
and inner planets, and the ratio of these eccentricities. 




5 10 15 

time/[10' P(1 AU)] 



Fig. 3. This figure shows the evolution of the planet semi- 
major axes and eccentricities for the run Rl shown in ta- 
ble |l|. 



In summary our simulation gives plasuble eccentric- 
ity values for the two planets, that can be understood in 
outline by use of a simplified analytic theory and are con- 
sistent with the current observations. 



4. Orbit integrations 

In addition to the analytic model presented in sec- 
tion (j^) and the hydrodynamic simulations presented in 
section (^, we have also performed three-body orbit in- 
tegrations using a fifth-order Runge-Kutta scheme (e.g. 
Press et al. 1993). 

The basic assumptions of the model are that the two 
planets exist within the inner cavity of a tidally truncated 
disc that lies exterior to the outer planet. Tidal interac- 
tion with this disc causes inwards migration of the outer 
planet, and also leads to eccentricity damping of the outer 
planet. It is further assumed that as the planets migrate 
inwards and approach their final semi-major axes, the disc 
disperses on some prescribed time scale tdisp- In our nu- 
merical calculations, a torque was applied to the outer- 
most planet such that it migrated inwards on a time scale 
oitjnig local orbital periods as defined in section(|^), and a 
damping force was applied in the radial direction to damp 
the eccentricity on a time scale of tc local orbital periods 
also as defined in section (||). 

These integrations used initial conditions correspond- 
ing to the more massive, outermost planet being located 
initially at 5 AU, with the lighter inner-most planet lo- 
cated initially at 2.5 AU. The planet masses adopted for 
the orbit integrations are the same as the minimum masses 
reported for the pl anets in the system around GJ876 by 
Marcy et al. ( |200l| ) (i.e. 1.87 Mj and 0.56 Mj). The stel- 
lar mass is taken to be 0.32 M©. Whilst these calculations 
provide only a crude approximation to the detailed physics 
of disc-companion interactions, their simplicity allows us 
to perform many calculations, covering a wide area of pa- 
rameter space, and also to run for much longer time scales 
than is possible for simulations of the type described in 
section(pl). 



A number of calculations have been performed to ex- 
amine the relationship between the final values of ei, 62, 
and their ratio ei/e2, to the various input parameters imigi 
tc, and tj^isp. The results of some of these calculations are 
presented in table and are discussed below. The unit of 
time used in the abscissa of the figures || to ^ is the orbital 
period for an object at 1 AU in orbit around a star with 
mass 0.32 Mq, and is denoted as P(l AU). 

4-.1. Dependence on migration and circularization times 

Equation ^ shows that the eccentricity of the outer 
planet, ei, depends on the ratio of tc/tmig- Here we present 
results of simulations that explore how the eccentricity ra- 
tio 62/61 depends on tc and tmig- 

Figure ^ shows the evolution of the semi-major axes 
and eccentricities for the run Rl, whose model parameters 
are described in table |^. This figure shows the inward mi- 
gration of the outer planet that subsequently locks to the 
inner planet as it reaches the 2:1 commensurability. The 
subsequent evolution is such that the two planets, now 
resonantly locked, migrate inwards. The eccentricities re- 
sult from the balance between eccentricity driving through 
the resonant interaction, and eccentricitiy damping due to 
the disc interaction. As the planets approach their final 
semi-major axes, the effects of migration and eccentricity 
damping are removed on a time scale of tdisp = 8 x 10'^ lo- 
cal orbits, causing them to cease migration at semi-major 
axes ai ~ 0.2 and a2 c± 0.126 which are values appropri- 
ate to GJ876. The subsequent evolution beyond a time 
of i ~ 0.45 P(l AU) in figure ^ occurs in the absence of 
disc effects, and suggests a long-term stability of the sys- 
tem given that it remains stable for 2 x 10^ orbits of the 
outer planet at its final semi-major axis. Figure ^ shows 
the evolution of the resonant angles (j) and tpi which librate 
about ^ = 1/) = in agreement with the results presented 
in section Bl for the full hydrodynamic simulations. 

Figure ^ shows the results from run R2 described in 
table |l], and illustrates the effect of reducing tmiq while 
keeping tc constant. As expected from equation Hq, the 
eccentricity of the outer planet increases from ei c± 0.1 to 
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<p v's t 



■f v's t 




0.5 1.0 1.5 

T;fre/[10' P(l AU)] 



0.5 1.0 1.5 

Time/[10' P(1 AU)] 



Fig. 4. This figure shows the evolution the resonant angles 
(j) and ip for the run Rl described in table |l|. These show 
libration around = "0 = 0, in agreement with the results 
from the hydrodynamic simulation shown in figure |l|. 



ei ~ 0.35 when tmig changes from 3.33 x 10^ to 3.33 x 10^ 
local orbital periods. However, the ratio 62/61 has also 
changed from 62/61 ~ 4 to 62/61 ~ 2. 

Calculation R3 illustrates the effect of keeping tmiq 
constant while changing tc- As expected from equation nS, 
reducing tc leads to a reduction in 61 , since the disc model 
damps the eccentricity more effectively. We also find that 
the ratio 62/61 again changes, it now being 62/61 ~ 6, as 
compared to 62/61 ~ 4 for run Rl. 

Equation ^ predicts that keeping the ratio tmig/tc 
constant, but changing both tmig and tc independently, 
should leave 61 unchanged. Calculation R4 indicates that 
this is what happens, and also shows that the ratio 62/61 
remains unchanged. 

Overall the results are entirely consistent with the 
analytic predictions presented in section ^ and with the 
hydrodynamical simulations presented in section ^. Fur- 
thermore, they indicate that the ratio 62/61 scales rather 
weakly with tmig/tc- We comment that from section ^ 
equation 10 we expect the eccentricity ratio to reach ~ 11 
as 61 — > 0. Long-term stabilty of two-planet systems that 
become locked due to disc-induced orbital migration is 
also indicated by our calculations. In particlular run Rl 
covers a time scale corresponding to 2 x 10*" orbits of the 
outer planet in its final configuration. We find that it is 
possible to arrange ei to match the observed value of the 
outer planet in the GJ876 system by fixing tmig and choos- 
ing tc appropriately, but it is difficult to then obtain a 
value for 62 that matches the reported value of 62 = 0.27. 



4-2. Dependence on Disc Dispersal Time Scale 

We have performed simulations to examine the effect of 
removing the disc on different time scales on the stability 
of the system. We find that the stability is largely unaf- 
fected by the rate at which the disc is removed. Figure || 
shows the evolution of the semi-major axes and eccentric- 
ities from a run in which the disc was removed on a time 
scale of tdisp = 8 X lO'^ local orbital periods. In this case 
the disc dispersal was switched on once the outer planet 



semi-major axis was ai < 0.5 AU. The final semi-major 
axis of the planet is determined by the radius at which the 
disc dispersal is initiated, and the time scale over which 
the disc is removed. Calculation R5 was similar to Rl 
except that the disc was removed instantaneously. The re- 
sults presented in table |l| show that this has little effect 
on the final outcome. A slight increase in the scatter of 
the temporal distribution of eccentricities was observed. 



Semlmojor Axis v's Time 



Eccentricity v's Time 




0.00 0.05 0.10 0.15 0.20 0.00 0.05 0.10 0.15 0.20 
linne/[10' P(l AU)] time/[10' P(1 AU)] 



Fig. 5. This figure shows the evolution of the planet semi- 
major axes and eccentricities for the run R3 shown in ta- 
ble |l|. 



5. Summary and Discussion 

In this paper we have considered two protoplanets gravita- 
tionally interacting with each other and a protoplanetary 
disc. The two planets orbit interior to a tidally maintained 
disc cavity while the disc interaction induces inward mi- 
gration. 

We have supposed that the previous evolution of the 
system results in the the planets getting into a configu- 
ration with an orbital separation just exceeding that re- 
quired for a 2 : 1 commensurability. This evolution is likely 
to have involved both accretion and migration. Given the 
set up we consider a natural evolution is towards both 
planets orbiting in a cavity outside of which orbits the 
protoplanetary disc. Tidal interaction results in the in- 
ner disc material either being expelled into the outer disc 
or accreting onto the central star. Subsequently the outer 
planet migrates towards the inner one as a result of inter- 
acting tidally with the exterior disc. 

When the migration is slow enough, we found that 
the outer protoplanet approached and became locked into 
a 2 : 1 commensurability with the inner one. This was 
maintained in subsequent evolution. We studied the na- 
ture of these interactions using a simple analytic model, 
hydrodynamic 2D simulations and longer time N body in- 
tegrations. These all gave consistent results. 

The magnitude of the stabilized eccentricities was 
found to be determined by the ratio of the migration rate 
to the circularization rate induced in the outer planet or- 
bit by the external disc. The eccentricity ratio 61/62 was 
found to vary with the magnitude of the ratio being more 
extreme for smaller eccentricities. 
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We have applied our results to the recently ciiscovcrcd 
resonant planets around GJ876. Simulation shows that a 
disc with parameters expected for protoplanetary discs 
causes trapping in the 2 : 1 commensurability when the 
planets orbit in an inner cavity and that eccentricities in 
the observed range may be obtained. In this case the or- 
bits were found to be aligned with both resonant angles 
librating about zero. The whole system then precessed in 
a retrograde sense. Finally when the disc is removed on 
a range of timescales the orbital configuration has been 
found to be stable for up to 2 x 10^ orbits subsequently. 
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